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The intensity of low frequency gravitational waves from black hole MACHO 
binaries is studied. First we estimate the gravitational wave background pro- 
l/"") ' duced by black hole MACHO binaries in the Milky Way halo as well as the 

cosmological gravitational wave background produced by the extragalactic 
Q\ ■ black hole MACHO binaries. It is found that the cosmological gravitational 

wave background due to black hole MACHO binaries is larger than the halo 

00 . 

\ background unless an extreme model of the halo is assumed, while it is smaller 

(— I than the background due to close white dwarf binaries at Ug W < 10" 2 - 5 Hz if 

Qh' the actual space density of white dwarfs is maximal. This cosmological back- 

er ' ground due to black hole MACHO binaries is well below the observational 

constraints from the pulsar timing, quasar proper motions and so on. We find 
that one year observation by LISA will be able to detect gravitational waves 
from at least several hundreds of nearby independent black hole MACHO bi- 
naries whose amplitudes exceed these backgrounds. This suggests that LISA 
j_j ■ will be able to pin down various properties of primordial black hole MACHOs 



together with the results of LIGO-VIRGO-TAMA-GEO network. Further- 
more, it may be possible to draw a map of the mass distribution of our halo, 
since LISA can determine the position and the distance to individual sources 
consisted of black hole MACHO binaries. Therefore, LISA may open a new 
field of the gravitational wave astronomy. 

PACS numbers: 98.80.-k, 04.30.-w, 97.60.Lf, 95.35.+d 
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I. INTRODUCTION 



The observations of gravitational microlensing toward the Large Magellanic Cloud (LMC) 
have revealed that a significant fraction of the Milky Way halo consists of solar or sub-solar- 
mass compact objects, which are called massive compact halo objects (MACHOs). The 
analysis of the first 2.1 years of photometry of 8.5 x 10 6 stars in the LMC by the MACHO 
Collaboration Jj]] suggests that the fraction 0.62^2 of the halo consists of MACHOs of 
mass 0.5^_q^Mq assuming the standard spherical flat rotation halo model. The preliminary 
analysis of four years of data suggests the existence of eight additional microlensing events 
with tdur ~ 90 days in the direction of the LMC 0. 

At present, we do not know what MACHOs are. This is especially because there are 
strong degeneracies in any microlensing measurements among the mass, the velocity and 
the distance to a lens object. There are several candidates proposed to explain MACHOs. 
However, all of them have some theoretical or observational difficulties. The inferred mass is 
just the mass of red dwarfs. However, observationally there is a tight constraint on the mass 
fraction of red dwarfs to the total mass of the halo [0-0]. Red dwarfs can contribute to, at 
most, a few percent of the total mass of the halo. The possibility that the MACHOs are neu- 
tron stars is ruled out by the observational constraints on the metal and helium abundance 
|J. If MACHOs are white dwarfs, the initial mass function (IMF) must have been peaked 
around ~ 2M Q when they formed [0-0, which IMF is completely different from the present 
day IMF. Furthermore, assuming the Salpeter IMF with an upper and a lower mass cutoff, 
the mass fraction of the white dwarfs in the halo should be less than 10 % from the number 
count of high z galaxies [l(| . The observation of the chemical yield also disfavors MACHOs 
being white dwarfs |TT],12|. Brown dwarfs are possible candidates for MACHOs [13 



since 



they are free from problems concerning metals and star counts. However, they require both 
a non-standard halo model and a mass function concentrated close to the hydrogen burning 
limit. Extrapolating the mass function of red dwarfs, we find that the contribution of brown 
dwarfs to the mass of the halo is less than a few percent 11 . In any case, extreme param- 



eters or models are needed as long as white dwarf/brown dwarf MACHOs are concerned, 
although future observations of the high velocity white dwarfs/brown dwarfs in our solar 
neighborhood might prove the existence of white dwarf /brown dwarf MACHOs (see also 

0)- 

Measurements of the optical depth in other lines of sight, including SMC and M31, are 
needed to confirm that MACHOs exist everywhere in the halo. One microlensing event 
toward SMC ||T6| , [iT|| is not sufficient to determine the optical depth toward SMC reliably. 
Since only the optical depth toward the LMC is available at present, in principle there 
is the possibility that MACHOs do not exist in other lines of sight. Any objects clustered 
somewhere between the LMC and the sun with the column density larger than 25M Q pc -2 |18| 



can explain the data. They include the possibilities; LMC-LMC self-lensing, the spheroid 
component, the thick disk, a dwarf galaxy, the tidal debris, warping and flaring of the 
Galactic disk [TP-^]. However, none of them do not convincingly explain the microlensing 
events toward the LMC. This means that it is worth considering alternative possibilities 
that MACHOs are not stellar objects but absolutely new objects such as black holes of mass 
~ 0.5M Q or boson stars with the boson mass of ~ 10~ 10 eV. 

In this paper, we consider primordial black holes to explain MACHOs. This possibility is 
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free from observational constraints at present. Theoretically it is still uncertain if black holes 
are formed in the early universe in spite of several arguments on the formation mechanisms 
of primordial black holes |23|-f25|j. Our standpoint here is not to discuss the details of the 



formation mechanism but to establish the observational signatures of primordial black holes 
if they exist. Firstly, electromagnetic radiation from gas accreting to a black hole MACHO 
(BHMACHO) is too dim to be observed unless the velocity of BHMACHOs is exceptionally 
small []26[ |. If primordial black holes are formed, most of them evolve into binaries through the 
three body interaction in the early universe at t ~ 10 _5 s p7| , |28[| . Some of these BHMACHO 
binaries coalesce at or before present epoch and it may be possible to detect the gravitational 
waves from such coalescing BHMACHO binaries. The event rate of coalescing BHMACHO 
binaries is estimated as 5 x 10~ 2 x 2 ±:L event s/yr /galaxy, which suggests that we can detect 
several events per year within 15 Mpc by LIGO-VIRGO-TAMA-GEO network [0,0. 

On the other hand, the other part of BHMACHO binaries still have large separations, 
and they emit gravitational waves at lower frequencies. These low frequency gravitational 
waves may be detected by the planned Laser Interferometer Space Antenna (LISA), which 
will be sensitive to gravitational waves in the frequency band of 10~ 4 — 10 _1 Hz. In this 
paper, we investigate low frequency gravitational waves emitted by BHMACHO binaries]] 
and we argue that LISA will be able to detect gravitational waves from BHMACHO binaries. 

The mass of MACHOs inferred from observations depends on the halo model (e.g. 



Ref. [0). Here we simply adopt M ~ 0.5M Q as a possible mass of MACHOs follow- 
ing the suggestion by MACHO Collaboration, and then investigate the dependence of the 
results on the mass. 

In Sec. II we review the formation scenario of BHMACHO binaries proposed in Ref. 



27| , |28|| . In Sec. Ill we establish a method to estimate the probability distribution function 
of the binary parameters as a function of the age after the binary formation, t. In Sec. IV we 
calculate the gravitational wave background produced by BHMACHO binaries in the Milky 
Way halo. In Sec. V we obtain the cosmological gravitational wave background produced 
by extragalactic BHMACHO binaries. We also compare the amplitude of the cosmological 
gravitational wave background produced by BHMACHO binaries with the observational 
constraints such as the pulsar timing and quasar proper motions. In Sec. VI, taking into 
account other noise sources in the LISA detection band, such as instrumental noises and the 
background produced by close white dwarf binaries (CWDBs), we summarize the spectral 
density of all noises in terms of the gravitational wave amplitude. In Sec. VII we estimate 
how many nearby BHMACHO binaries are expected to be observed by LISA. Section VIII 
is devoted to summary and discussion. 



II. FORMATION SCENARIO OF BHMACHO BINARIES 

We briefly review the formation scenario of BHMACHO binaries given in Ref. p7|j28 



to introduce our notations. For simplicity, we assume that black holes dominate the dark 
matter, i.e., Q = Qbh, where Qbh is the density parameter of BHMACHOs at present. 



*W. A. Hiscock has also independently carried out a similar study using a simpler model. 
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Although black holes dominate the dark matter at present, their fraction of the total energy 
density is small ~ 10~ 9 at their formation from Eq. (|2.1|) , since the ratio of the radiation 
density to the black hole is in proportion to the scale factor. They are harmless to the 
nucleosynthesis when their fraction is only ~ 10~ 6 . BHMACHOs behave like cold dark 
matter for the large scale structure formation, and the large primordial density perturbations 
that are needed for BHMACHO formation are compatible with the observed upper limit of 
CMB spectral distortions [27]. Furthermore, we assume that all black holes have the same 



mass Mbh- We normalize the scale factor such that R — 1 at the time of matter-radiation 
equality. 

Primordial black holes are formed when the horizon scale is comparable to the 
Schwarzschild radius of the black holes [Q. The scale factor at the time of the black 
hole formation is given by 

Ri 



GM BH _ ,„ jMn,,^ 



where H eq is the Hubble parameter at the time of matter-radiation equality, i.e., cH~ q = 

J '3c 2 / 8tiG p eq = 1.2 x 10 21 (Qh 2 )~ 2 cm. The mean separation of black holes at the time of 
matter-radiation equality is given by 

x = (M BH /p eq )^ 3 = 1.2 x lQ l \M BH /M Q fl\nh 2 Y^ [cm]. (2.2) 

Consider a pair of black holes with the same mass M B h and a co moving separation x < x. 
We define an averaged energy density of this pair of black holes by dividing their masses by 
the volume of a sphere whose comoving radius equals to x as p BB = p eq x 3 / (x 3 R 3 ) . Then, 
p BH becomes larger than the mean radiation energy density, p B = p eq /R A , when 

x x 3 



R>R m = . (2.3) 

After R = R m , the motion of the pair is decoupled from the cosmic expansion, and the 
pair forms a binary. Unless x is exceptionally small, this newly formed binary is kept from 
colliding with each other by the effect of tidal force from neighboring black holes, which 
gives the binary sufficiently large angular momentum. 

The initial value of the semimajor axis a will be proportional to xR m , and hence 

a = axR m = a—, (2.4) 
x A 

where a is a constant of order unity. To estimate the tidal torque, we assume that the tidal 
force is dominated by the black hole that is nearest to the binary. We use y to represent the 
comoving separation of the nearest neighboring black hole from the center of mass of the 
binary. Then, the semiminor axis, b, will be proportional to (tidal force) x (free fall time) 2 , 
and hence it is given by 

M BH xR m {xR m f n(A 3 ( . 

{yR m y M BH \y I 



4 



where (3 is a constant of order unity. Then, the binary's eccentricity, e, is evaluated by 



\ 



1-p 2 




(2.6) 



In Ref. [p8[l , it is verified that the analytic estimates given in Eqs. (|2.4j ) and (|2.5| ) are good 
approximations and the numerical coefficients, a and /?, are actually of order unity. In this 
paper, we adopt a = 0.5 and (3 = 0.7. Q 

If we assume that black holes are distributed randomly, the probability distribution func- 
tion, P(x,y), for the initial comoving separation of the binary, x, and the initial comoving 
separation of the nearest neighboring black hole from the center of mass of the binary, y, is 
given by 



P(x, y)dxdy 



x 6 



dxdy, 



(2.7) 



where < x < y < oo. This probability distribution function is normalized to satisfy 
/ °° dx f£° dyP(x,y) = 1. Changing the variables x and y in Eq. ( |2.7|) to a and e with 
Eqs. (|2.4j) and (|2.6| ), we obtain the probability distribution function for the eccentricity and 
the semimajor axis of binaries as 



f(a, e)dade 







I.V2, 



4 (ax) 3 / 2 (1 - e 2 ) 3 / 2 



exp 







(1 _ e 2)l/2 



x 3/4 



dade, 



(2.8) 



where a/1 - /3 2 < e < 1 and < a < oo. Note that /(a, e) = for < e < VI - P' 2 - 

We consider here BHMACHO binaries whose semimajor axis is less than the mean sepa- 
ration x. Their coalescence time due to the emission of gravitational waves is approximately 
given by Pi 



t = tn 



fl 



'l-e 2 )i. 



a = 2.0 x 10 



ii 



M. 



BH 



cm 



(2.9) 



where t = 10 yr and ao is the semimajor axis of a binary in a circular orbit which coalesces 
in t . Using Eqs. ( |2.4| ) and ( |2.6| ), Eq. ( ^.9| ) can be written in terms of x and y as 



t 



-21 



0, _ 
a 



(2.10) 



Integrating Eq. ( p.7|) for a given t with the aid of Eq. Q2.10 ), we obtain the probability distri- 
bution function for the coalescence time ft{t)- We should take the range of the integration 
to satisfy < x < x and x < y < oo. The first condition (x < x) is necessary for the pair 



2 In Ref. pifl , a = 0.4 and (3 = 0.8 are adopted. However, these values are obtained by a least 
squares fitting without fixing the power indices of a and b. When we use the power indices of the 
analytical estimates, i.e., 3 for a and 4 for 6, it is better to adopt a = 0.5 and (3 = 0.7. Note that 
the results of this paper are not affected so much by the detailed choice of these constants. 
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to be a binary. The second condition turns out to be (t/f) 1 / 16 ^ < y < (t/t) l / 21 x for a fixed 
t. Performing the integration, we obtain 



3 / A 3/37 



58 ft^ 



37' \t 



r 



58 m- 1 / 7 



37' \t 



dt 
T 



3 ft\^ 37 /58\ dt 
\37J J'' 



37 W 



where T(x, a) is the incomplete gamma function defined by 

T(x, a) = 



(2.11) 



(2.12) 



and T(58/37) = 0.890- • •. The second equality in Eq. (|2.11| ) is valid when t/t <C 1. This 
is automatically satisfied as long as we consider the case t ~ to with typical values of 
parameters, for which to/i 1 is satisfied. 



III. LATE TIME DISTRIBUTION OF THE BINARY PARAMETERS 

A. numerical approach 

Since a BHMACHO binary radiates gravitational waves, the binary parameters change as 
a function of time due to radiation reaction. The probability distribution function, f(a,e), 
in Eq. Q2.8Q is also a function of time. In this section, we develop a method to compute this 
time dependent distribution function, /(a, e;ii), where t x is the time after the formation of 
binaries. We refer to the initial values of a and e as and e^, respectively. 

Firstly, we represent a, and as a function of the present values of binary orbital 
parameters, a and e, and the component masses of the binary, M BH i and M BH2 . The 
formula that relates to e is given in Ref. by 



19Bh pdee 29 / 19 [l + (121/304)e 



12c^ 



2l 1181/2299 



[1 - e 2 ) 3 / 2 



(3.1) 



where 



B 



64 G 3 M BH1 M BH2 (M BH1 + M BH2 ) 



(3.2) 



and 



-12/19/. 2^ 



121 „■ 

1 + m ei 



-870/2299 



ae 



-12/19 (1 _ e2l 



, 121 2 

1 H e 2 

304 



-870/2299 



(3.3) 



Using the above relations, we can determine cij and as a function of a, e and t% numerically. 
We can also obtain a* and e^, by numerically integrating the evolution equations 



da 
~dt 



B 



a 3(! _ e 2)7/2 



, 73 2 37 4 
1 + — e 2 + — e 4 



24 



96 



de 



19 



Be 



121 \ 

12a 4 (l-e 2 ) 5 / 2 V 1 + 304 6 J ' (3 ' 4) 
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from the present time backward to the initial time. 

Secondly, we need to calculate the Jacobian to relate the distribution of a and e with 
that of a, and e,. For this purpose, we take the total differentiation of Eqs. fl3.1| ) and ( p.3|) 
to obtain 

da 12 1 + (73/24)e 2 + (37/96)e 4 de da { 12 1 + (73/24)ef + (37/96)ef de; 



a 19 (1 - e 2 ) [1 + (121/304)e 2 ] e m 19 (1 - ef) [1 + (121/304)e 2 ] d 



da 121 + (73/24)e 2 + (37/96)e 4 de 
~~a ~ 19 (1 — e 2 ) [1 + (121/304)e 2 ] ~~e 



e 2 9 /i9 [i + (i21/304)e 2 ] 1181/2299 de 
+ (1 - e 2 ) 3 / 2 J 

e 29/19 [l + (121/304)e 2 ] 1181/2299 de, 
(1 - e 2 ) 3 / 2 J 



(3.5) 



Then, by taking the wedge product of the above two equations, we easily find the Jacobian 
relation, 

da e 29 / 19 [l + (121/304)e 2 ] 1181/22 " ^ _ da, e 29/19 [1 + (121/304)e 2 ] 1181/2299 

(1-e 2 ) 3 / 2 6 " a, (1-e 2 ) 3 / 2 ^ (3 ' 6) 

Since the distribution function f(a, e; t\) at a time ti is related to the initial distribution 
function /j(a;,ej) by 

/(a, e; ti) da de = /j(aj, e,) da, de i; (3.7) 

we find 



e \ 29/19 /i_ e 2\ 3/2 [l + (121/304)e 2 



1 + (121/304)e 2 



1181/2299 



fi(a, h ei), (3.8) 



where ^(a^e,) is given by Eq. fl2~lf) . 

The Keplerian orbital frequency of a binary, u p , is related to the length of the semimajor 
axis, a, by 



9 G(M Bm + M BH2 ) 

p = V a 1 ' ( } 

where we add the subscript, p, to the orbital frequency, u, for later convenience. Hence, the 
distribution function for a and e is transformed into that for u p and e by 



fv,e,t{^^t x ) = /(a,e;ti) 



da 



dz/p 



^f{a^h). (3.10) 



In Fig. 1, we plot f u , e ,t(v p , e; ti) for ti = to = 10 10 yr as a function of e for several orbital 
frequencies. Here we adopt M BH = O.5M and Qh 2 = 0.1. As we can see in Fig. 1, binaries 
with a high orbital frequency (v p <^ 10 -3 Hz) are almost in circular orbits at present, while 
binaries with a low orbital frequency [y v % 10 -3 Hz) keep their eccentricity large even now. 
The transition frequency from an eccentric orbit to a circular one depends on M B Hi Qh 2 
and £i. Note that almost all binaries formed through the mechanism presented in Section. |I| 
are eccentric when they are formed, as we can see from Eq. (|2[ 
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B. approximate formula 



An approximate formula for the distribution function of u p can be obtained in the case 
that the eccentricity has already decayed at a time t\. Under the assumption e < 1, the 
remaining time before a binary with the semimajor axis, a, coalesces is given by [|3T 



425 a 4 

At = —t — (3.11) 
768 a 

where a is given in Eq. (|2.9| ) and t = 10 10 yr. On the other hand, the distribution function 
of coalescence time is given by Eq. fl2.11| ), which can be interpreted as the distribution 



function of At by replacing t with ti + At. Then, by using Eqs. ( |3.9|) and (|3.11|) , we obtain 
an approximate formula for the distribution function of the orbital frequency as 

425 ^i\ 3/37 t ( a {58\ dv v 



where we assume t\ 3> At. This expression corresponds to the quantity that is obtained by 
integrating the distribution function f v ^A v pi e ! *i) * n Eq- (|3.10|) over e. Since the semimajor 
axis, a, is proportional to v' 2 ^ in Eq. (^.9|), we find 

f v ,t{ypM)<xv- n /*. (3.13) 

On the other hand, we have oc from Eq. (fT^), a oc 

M^h from Eq. and 

t oc M^ /3 (ft/i 2 ) 16 / 3 from Eqs. (P|), Q) and ( g3D| ). Then, for a fixed orbital frequency, 



we find there is a relation 

M^h)ocM B ]J 0/lu (nh 2 r/ 37 . (3.14) 

Since the mass of MACHOs has not been fully determined by the microlensing observation, 
we will consider the dependence of the results on the mass using this relation. 



IV. GRAVITATIONAL WAVE BACKGROUND PRODUCED BY BHMACHO 
BINARIES IN THE MILKY WAY HALO 

In this section, we investigate gravitational wave background produced by BHMACHO 
binaries in the Milky Way halo. In our model, there are so many BHMACHO binaries in 
our galaxy that they themselves become stochastic background sources, and may limit the 
confusion noise level for individual sources. 



A. gravitational waves from two point masses in a Keplerian orbit 

A circular orbit binary emits gravitational waves only at a frequency v gw = 1v v . However, 
a binary with eccentricity e ^ emits gravitational waves at frequencies v gw = pv p (p = 
1,2, ■ ■ •) [|32],[33|. As we have seen at the end of Section |III Aj , almost all binaries formed 



through the mechanism explained in Section || are highly eccentric when they are formed, 
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and orbits of binaries with a low orbital frequency still remain eccentric even at present. 
Therefore, the contribution from high harmonics may be important. 

Gravitational waves at frequency v gw are emitted by binaries whose present Keplerian 
frequency is 

Vp = Vgw/P, (p=l,2,---). (4.1) 

The gravitational wave luminosity from a binary which emits gravitational waves at fre- 
quency v gw as the p-th harmonic is given by |32] 



'gwi W ■'-'Wgw 

where 



LV>( Vgwt e) = L vl°J 3 p- 10 / 3 g(p, e), (4.2) 



L = 1.351 x 10- p^^AM. gl + M^ 2) y ^M BH1 + M BH2 y ^ Hz _ 10/3] 

(4.3) 

and Mbhi and M B H2 are component masses of the binary. The function g(p, e) is given by 



p4 

9M = -- 



2 1 
J p -2(pe) - 2eJ p -i{pe) + -Jpipe) + 2eJ p+1 (pe) - J p+2 (pe) 



+ (1 - e 2 ) [J p _ 2 (pe) - 2J p (pe) + J p+2 (pe)Y + — [J p (pe)Yj, (4.4) 
where J n (x) is the Bessel function. We note that g(p, 0) = for p ^ 2. 

B. optical depth toward the LMC 

Let us consider the density profile of the halo. Based upon the number and duration 
of gravitational microlensing events fl], the optical depth of MACHOs toward the LMC is 
estimated as 

r LMC ~ 2 x 1(T 7 . (4.5) 

While, the optical depth can also be expressed in terms of the number density of MACHOs. 
Let us assume that the distribution of the number density of MACHOs is spherical symmetric 
as 

[1 + {r 2 /Dl)\ x 

where r is the distance from the center of our galaxy, n s is the density of MACHOs at the 
galactic center, and D a is the core radius of the distribution |34j]. Then the optical depth is 
given by 

t lmc = ^GM bh r». / _ x\ {r)dX) (4.7) 



c 2 Jo V D 
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where x is the distance from the earth toward a lens, and D s is the distance to the LMC 
from the earth. Here we adopt D s = 50 kpc. In general, r can be expressed by x and the 
directional cosine, cos 9, as 



D 2 



2D x cos 9 + x 2 



(4.8) 



where D is the distance from the galactic center to the earth, and we adopt D = 8.5 kpc 
here. When we consider the LMC direction, we adopt cos 9 = rj := 0.153. If we fix A and 
D a , Eqs. (|4.5|) and (Of) determine n s as 



C 2 T LMC 

AttGMbhDI 



Da 
D„ 



dy 



u+(§ 



')}' 



(4.9) 



where the quantity in the square bracket is dimensionless. 



C. flux and strain amplitude of the halo background 



The contribution to the average flux from the p-th harmonic per unit frequency at v gw 
is given by 



" X ; n(r) I deL ip) (u gw ,e)f u ^t(iy p ,e;to)-^-[eTg cm 2 s 1 Hz x l 



And 2 



dv, 



(4.10) 



yw 



where d = (Dq — 2rD cos 9 + r 2 ) 1 / 2 is the distance between a source and the earth. The 
probability distribution function f u ,e,t{v p , e; to) is given by Eq. ( |3.10|) and the luminosity from 



the p-th. harmonic, L^(u gw ,e), is given by Eq. (O). Defining 

d 3 x 



I :-- 



And 2 



n(r), 



the dependence of the flux on the halo model is factorized as 

Fjf\v gw ) = IL u l g °J 3 p- 13 / 3 £ def Ute>t {v p , e; t )g(p, e) [erg 
Using Eq. (|J) with Eq. (p|), we find 



cm" 2 s" 1 Hz" 1 ! 



(4.11) 



(4.12) 



C 2 T LMC 



16nGM BH D c 



Q 2 T LMC 



16nGM BH D c 



ydy 



{ (v-a)V (1+ ^. 



Ds 



dy 



y{ l -ir s y 



2 U 



where we defined the halo shape factor, /, by the last equality. Note that / is a dimensionless 
constant which depends only on the shape parameters of the halo (A and D a ). In Fig. 2, we 
plot J as a function of the core radius, D a , for several values of A. As we can see in Fig. 2, 
/ is about 7 when the halo of our galaxy is isothermal, i.e., when A = 1. If the halo is more 



(4.13) 
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concentrated to the galactic center, i.e., A is larger and D a is smaller, the gravitational wave 
flux becomes larger. 

Collecting the contributions from all harmonics, we obtain the total gravitational wave 
flux per unit frequency at v gw as 



p=l 



oo ,i 

= IL v 1 g °J 3 J2P~ 13/3 defu,e,t(.Vp,e;to)g(p,e) [erg cm" 2 s" 1 Hz" 1 ]. (4.14) 

When we evaluate the above expression numerically, we took the summation over p up to 
p = 1000. In Appendix |A], we show that the error caused by neglecting higher p harmonics 
is less than a few percent. 

For comparison, it is convenient to translate the gravitational wave flux into the spectral 
density of the gravitational wave strain amplitude, h u , which is calculated by 



h h u al °= |^jL— = 5.615 x 10- 20 ^! [Hz" 1 / 2 ]. (4.15) 



" 3 / 2 a/7T U gw U gw 



In Table I, we list the spectral density h„ /v I at several frequencies for several BHMACHO 

t2 M^+^fV,of hhalo , 



masses, M BH , and density parameters, Qh . Note that h„ / V / does not depend on the halo 
shape factor, /. In Fig. 3, the solid line is the spectral density h^ alo /\/l for Mbh = O.5M 
and Qh 2 = 0.1. We also plot the spectral density obtained by taking the summation up 
to p — 10 in Eq. ( 4.14 ) by the long dashed line. We see that the contribution from higher 
harmonics takes maximum at about 10 -3 Hz. Even there, the contribution is not so large, 
and the results coincide with each other within 25%. 



D. approximate relations for the halo background 

When the condition e <C 1 is satisfied, a good approximation for the halo background can 
be obtained analytically. Since only the second harmonic contributes to the gravitational 
wave flux in this case, the total gravitational wave flux in Eq. ( |4.14| ) can be approximated 
by 

F v {v gw ) ~ /Lo^f 2" 13/3 /,> P ; *o) : e « 1, (4.16) 

where fv,t{v p \ t) is given by Eq. ( |3.12| ). Since a relation F v {y gw ) oc v g ^ 3 can be derived from 
Eqs. (|3.13p and (|4.16|) , the spectral density h^ al ° in Eq. ( |4.15|) satisfies 

K oc v~ll % : e < 1. (4.17) 

On the other hand, for a fixed gravitational wave frequency, we have I oc M^jj- from 
Eq. ( P3D , L oc M^ 3 from Eq. Q), and f v>t (v p ; h) oc M^ 70/m (fi/i 2 ) 16 / 37 from Eq. {^J§. 
Then, by using Eqs. (|4.15|) and ( |4. 16|) , we find 

h h „ al °otM 8 B % 222 (nh 2 ) 8 / 37 :e«l. (4.18) 
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As we can see in Table I, these approximate relations hold for v gw > 10 _1 Hz, at which 
almost all binaries that contribute to the halo background are in circular orbits. 

In Fig. 3, we plot the approximate spectral density by the dotted line, replacing F v 
in Eq. ( }4.15|) with the approximate one in Eq. ( f4.16|) . The dotted line coincides with the 
solid line quite well at high frequencies, i.e., v gw > 1CT 2 Hz. This means that almost all 
binaries that contribute to the gravitational waves at high frequencies are in circular orbits. 
Conversely, the fact that the spectral density h^, al ° deviates from the approximate value at 
low frequencies means that binaries with large eccentricity e contribute the gravitational 
wave flux at low frequencies. 



V. COSMOLOGICAL GRAVITATIONAL WAVE BACKGROUND 

In the previous section, we studied the gravitational wave background produced by BH- 
MACHO binaries in the Milky Way halo. In this section, we compute the gravitational 
wave background produced by extragalactic sources, i.e., the cosmological gravitational wave 
background, which turns out to dominate the background radiation. 

The gravitational wave flux is proportional to the factor, I = j d 3 xn(r) /And 2 ~ n(D )r. 
Strictly speaking, this expression needs some caution because we need to take into account 
the expansion of the universe when we estimate the cosmological background. However, as 
a rough estimate, we can use the estimate I ~ nrh orizon for the cosmological background, 
where n is the mean number density of BHMACHOs and rh orizon is the horizon radius. On 
the other hand, since the mean number density of BHMACHOs in the halo is enhanced 
by the factor (r galaxy /r ha i ) 3 , we have I ~ n(r gaia:r 2 / /r /iaio ) 3 r /iaio for the halo background, 
where r ga i axy is the mean separation between galaxies and rtaio is the radius of the halo. 
Then, the ratio of the cosmological background to the halo background is estimated as 
{r horizon /r ga iaxy){r halo /r g aiaxy) 2 ~ (3000Mpc/lMpc) (0.05Mpc/lMpc) 2 ~ 7.5. This suggests 
that the cosmological background can be the same order of the halo background and it 
cannot be neglected. 



A. cosmological model 

We assume that the cosmological distribution of BHMACHOs is homogeneous and 
isotropic. The Hubble equation is 

u2 [R\ 2 8ttG Kc 2 

:= \R ) =^ p ~^ (5 ' 1} 

where we assume that the cosmological constant is zero. For simplicity, we assume that the 
energy density is determined by BHMACHOs and the radiation. The energy densities of 
BHMACHOs and the radiation at present are denoted by pbh and pr, respectively. We refer 
to the present density parameter of each component as Om = Pbh/Pc and Qr = Pr/ p c , 
where p c = 3Hq/8tiG is the critical density and H = lOO/i km/s/Mpc. Since the energy 
density of BHMACHOs behaves like that of dust, the total energy density p is given by 

Q ~ = + T&Vr, tt = Q BH + Q R , (5.2) 
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where quantities with the subscript, 0, represent the present values. Using Eq. (|5.1| ), the 
relation between the scale factor and the cosmological time, t, is obtained as 



H dt 



(*)"U) 



n* + £W (£) + (l - !)) ) : 



1/2 ' 



(5.3) 



where we use Eq. (|5.2|) and the relation, 



Kc 2 



(5.4) 



Since the present temperature of the cosmic microwave background is about 2.75 K, the 
energy density of the radiation is given by = 4.81 x 10 _34 g cm -3 . Recalling that the 
scale factor is normalized such that R = 1 at the time of matter-radiation equality, we see 
that #0 is given by 1/R = Vt R /Vt BH = 2.56 x 1(T 4 (fWi 2 ) -1 ~ 2.56 x 1(T 4 {Vth 2 )' 1 . 



B. energy density and strain amplitude of the cosmological gravitational wave 

background 

When we consider the cosmological situation, we have to take the redshift into account. 
Let us consider gravitational waves observed at a frequency v gw . If they are emitted by 
a binary at a redshift 1 + z = Ro/R as the p-th harmonic, the frequency of them at the 
moment of emission is h> gw Ro/R, and the Keplerian orbital frequency of the binary is given 
by 

MR) = • ( 5 - 5 ) 

R p 

The gravitational wave luminosity emitted by this binary in the p-th harmonic is evaluated 

by@ 

LiP) (f^' 6 ) = Lo (t) 1073 ^ 10 ^^^), (5.6) 

where L is given in Eq. ( |4.3| ). Then, the gravitational wave luminosity per unit logarithmic 
frequency at v gw per BHMACHO is given by 

Caw{vm,,\t) = V / deL (p> ( ^-u m „,e] f,,. R /(z/„, e; tW, 



9W dv 



GW {Vgw] t) = J de LiP) (jJ[ U 9 w > e ) f^^P' e ' f ) 1 

/ d \ 13/3 oo „ 

= LwfJ* {-£) J2P~ 13/3 J deU, t (v P , e; t)g(p, e) [erg a" 1 ]. (5.7) 

Of course, the binary systems that coalesced at some time t < t are included, although 
the binary systems had not coalesced when they emitted gravitational waves observed at a 
frequency v gw by us. 

In order to characterize the spectrum of stochastic gravitational waves, we often use 
^Gw{ v gw)i the ratio of the energy density of the gravitational waves per unit logarithmic 
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frequency to the critical density f35f . To evaluate flawi^gw), we use the fact that the 
gravitational wave energy existing in a unit comoving volume is equal to the red-shifted 
total gravitational wave energy that is emitted per unit comoving volume until now. The 
gravitational wave energy per unit logarithmic frequency in a unit comoving volume is given 
by p c c 2 flGw( l/ gw)Ro- On the other hand, the latter quantity can be evaluated as follows. The 
comoving number density of BHMACHO is given by (1/x) 3 , where x is the mean separation 
at the time of matter-radiation equality given in Eq. fl2.2| ). Then, the red-shifted total 
gravitational wave energy that is emitted per unit comoving volume per unit logarithmic 
frequency during t ~ t + dt is given by Cqw {y gw ;t) (l/x) 3 dt. Integrating this expression 
multiplied by the red-shift factor (R/R ), we obtain the quantity that is to be equated with 
p c c 2 &gw (ygw ) Rq ■ 



From the above consideration, we obtain 
p c c 2 Jt f \RoJ \xR , 



tt G w(vgw) = — I dt ( — J ( — — J C GW (u gw ; t) 

' R : \ ° ( l 



^ | R ^ C G W (Vg W ;t(R)) \R ) [xR; ^ ^ 



R m /Ro \RqJ PcC 2 H 



n* + n OT (;j£) + (i-n) (f 



The relation between t and R is given by Eq. ( |5.3| ) and the binary formation time, tf, is the 
time at which R equals to R m given in Eq. ( |2.3|) . Hereafter, we set h = 0.8, when quantities 
with p c c 2 or Hq are evaluated numerically. In the case of h = 0.8 and Qh 2 = 0.1, the present 
age of the universe, to, is about 1.06 x 10 10 yr. 

In Fig. 4, we plot dQGw(v gw )/d()nR), which is the integrand of the last line of Eq. Q5.8D , 
as a function of the scale factor R/Ro by the solid lines for Mbh = O.5M and flBiih 2 = 0.1 
at several frequencies. When we evaluate C-Gw(y gw ]t) given in Eq. ( |5.7|) numerically, we 
perform the summation up to p = 10. The computational time does not allow us to sum up 
higher harmonics. Since the power index of dflGw(v gw ) / d(\ia R) with respect to R is larger 
than 0, the contribution to the energy density from binaries at higher redshift is smaller. 
Therefore, we perform the integration in Eq. (|5.8|) from R/Rq = 1/Rq (not R m /Ro) to 
R/Ro — 1 to evaluate the energy density of the cosmological background numerically. f\ 

In Fig. 5, we plot the density parameter of the cosmological gravitational wave back- 
ground per unit logarithmic frequency, ^Icwi^gw), for Mbh = 0.5M Q and flh 2 = 0.1. We 
again take the summation until p = 10 in Eq. ( |5.7| ). However, as we have seen at the end 
of Section. [IV C, the contribution from higher harmonics is not so large. For example, we 



calculated the density parameter of the cosmological gravitational wave background per unit 
logarithmic frequency at v gw = 10 -7 Hz for Mbh = 0.5M Q and Qh 2 = 0.1. When the sum- 
mation is taken until p = 10, we have flcwi^gw) = 2.7 x 10~ 20 . While, when the summation 
is taken until p = 1000, we have VLQ W {y gw ) = 3.2 x 10~ 20 . 



3 Note that, when the redshift is too high, the orbital frequency in Eq. Q5.5| ) may become larger than 
~ 1000 Hz. Since such binaries are not allowed to exist for black hole binaries with Mbh ~ Mq, 
the summation over them should be eliminated. However, such a complication is unnecessary as 
long as we consider v gm < 10 _1 Hz, for which at R > 1 the orbital frequency is less than ~ 1000 
Hz. 
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The spectral density of the cosmological gravitational wave strain amplitude is given by 



5.615 x Kr 20 -^ 3/2 [Hz~ 1/2 ]. (5.9) 

Vgw 

Note that the flux from the cosmological BHMACHO binaries per unit frequency at u gw 
is p c c 3 Qcw (^gw)/^gw, which corresponds to F v {y gw ) in Eq. ( |4.14|) for the halo background. 
In Table II, we list the spectral density, h c ° s , at several frequencies for several BHMACHO 



masses, Mbh, and the density parameters, flh . We present the figure of h c ° s in Sec. VI 



C. approximate relations for the cosmological background 

When the condition e <C 1 is satisfied, Eq. (|5.7|) can be approximated by 

, f? \ 13 / 3 

13/3 I £^0 \ 9 -13/3 



C GW (v gw ; t) ~ L u L g 6 J a ( ) 2- u ' 6 Ut{v p - 1) : e < 1, (5.10) 



where f v ,t{v p ] t) is given in Eq. fl3.12p . When Qbh is of order unity, R is approximately pro- 
portional to t 2 / 3 during the matter dominant era. For a fixed gravitational wave frequency 
u gw , we can extract independence from Eq. G3.12|) as f u ,R(v p ; R) oc t~ 34//37 z/ p (i?)~ 11//3 oc 



i? 254 / 111 , where we used the relation a oc v p (R)~ 2 l 3 oc R 2 ^ 3 that follows from Eqs. Q3.9| ) 
and (|5.5|). Hence, from Eq. ( |5.10| ) we can see that Ccwi^gw] t) oc R- 227 / 111 [ s approxi- 
mately satisfied. Since [fin + VL BH (R/ 'R ) + (1 — Vt)(R/ Rq) 2 }' 1 / 2 can be approximated by 
^l B ]^ 2 (R/ Rq)' 1 / 2 during the matter dominant era, we find that the integrand in Eq. ( |5.8|) 
is approximately proportional to R- 121 / 222 . i n Fig. 4, we also plot the approximation of 
dVLcwi^gw) I 'd(\nR) by the dotted lines, replacing the luminosity Caw in Eq. (|5/7|) with the 
approximation given in Eq. ( |5.10| ). The solid lines coincide with the dotted ones well for 
high frequencies. At high redshift almost all binaries that contribute to the gravitational 
wave flux at all relevant frequencies are in circular orbits, as is seen from the agreement of 
the solid and the dotted lines. 

When the approximation (|5.10| ) is valid, with the aid of Eq. (|3.13| ), it implies 



£Gw{vgw]t) oc v 2 J 3 . Thus we find Vt GW (v gw ) oc v 2 J 3 , and hence 

h™ oc v-V* : e « 1. (5.11) 
This dependence on v gw is the same as that obtained in the halo background case in 



Eq. ( 4.17 ). For fixed v gw and t, we have approximate relations (1/x) 3 oc M BH from 
Eq. (|3), L oc M B % 3 from Eq. and f v , R (v p ;R) oc M~^ 70/m from Eq. (^J§. There- 



fore, Qcwivgw) oc M B % 1U is derived from Eqs. (|5.8|) and ( |5.10|) , and hence we find there is 
an approximate relation 

CocMff 22 :e<l. (5.12) 

As we can see in Table. II, this relation is satisfied for v gw > 10 _1 Hz, at which almost 
all binaries that contribute to the cosmological background are in circular orbits. The 
dependence of the spectral density, h c ° s , on the density parameter Qh 2 cannot be derived so 
easily in an analytic method. 
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D. observational limit 



There are several observational constraints on the cosmological gravitational wave back- 
ground. These constraints come from the millisecond pulsar timing, quasar proper motions, 
the standard model of big bang nucleosynthesis, and so on. In this section, we study whether 
or not the cosmological gravitational wave background produced by BHMACHO binaries 
satisfies these constraints. Note that the frequency range considered in this section is much 
less than that of LISA. 

At the time of the nucleosynthesis, the effective number of neutrino species is restricted 
to be less than three, to be consistent with the observations of He, D, and Li 1311. If the 



effective number of neutrino species were larger than three, the expansion rate at the time 
of the nucleosynthesis would become larger and the predicted abundance of He, D and Li 
would differ. Since the gravitational wave background contributes to the energy density 
in the same way as neutrinos, the constraint on the effective number of neutrino species 
restricts the energy density of the cosmological gravitational wave background within that 
of one massless degree of freedom in thermal equilibrium. If the cosmological gravitational 
wave background already exists at the time of the big bang nucleosynthesis, ^Icwi^gw) is 
restricted approximately by 

d{\nv gw )tl GW {v gw )h 2 <HT 5 . (5.13) 

Here, we should note that gravitational wave background produced after the nucleosynthesis 
is not restricted by Eq. ( |5.13| ) at all. 



Observations of stable millisecond pulsars provide a constraint on the energy density 
of the gravitational waves |37 , |38| . Gravitational waves deform the metric perturbing the 



observed pulsar frequency. Therefore, the fluctuations in pulse arrival times are used to 
constrain the gravitational waves. Pulsar timing residuals are most sensitive to gravitational 
waves at frequencies around 1/T, where T is the total data span. Observations of PSR B 
1855+09 yield an upper limit, 

&Gw(vgw)h 2 < 1.0 x 10~ 8 (95% confidence), (5.14) 

at v gw ~ 4.4 x 10~ 9 Hz. 

After taking into account the energy carried away by gravitational waves and the effects 
of Galactocentric acceleration, orbital periods of binary pulsars are sensitive to gravitational 
waves with periods as great as the light-travel time from the pulsar. Current limits from 
such data are [j38| 

Vt GW {v gw )h 2 < 0.04 (95% confidence), (5.15) 

in the range 10~ n Hz < v gw < 10~ 9 Hz, and 

VL GW {u gw )h 2 < 0.5 (95% confidence), (5.16) 

in the range 10 _12 Hz < v gw < 10 _11 Hz. 

Background gravitational waves randomly scatter photons from distant quasars on its 
way to the earth. This may cause quasar proper motions. Using measurements of the proper 
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motions of quasars, upper limits on the energy density of gravitational waves with periods 
longer than the time span of observations is obtained [[39 



as 



tt GW {vgw)h 2 < 0.11 (95% confidence), (5.17) 

in the range u gw < 2 x 10~ 9 Hz. 

However, these limits are not enough to reject the existence of BHMACHOs. As we 
can see in Fig. 5, the cosmological gravitational wave background produced by BHMACHO 
binaries, ^Gwi^gw), is less than these limits at 10~ 9 Hz < v gw < 10 _1 Hz. At higher 
frequencies v gw > 1CT 1 Hz, ^lawi^gw) continues to grow as Vtcw^gw) oc v 1w till v gm ~ 10 3 
Hz. However, the constraint in Eq. (|5.13|) is still satisfied even at frequencies v gw ~ 10 3 
Hz. Furthermore, since almost all background at these frequencies is produced after the 
nucleosysthesis, Eq. ( |5.13 ) is not a real constraint for the existence of BHMACHOs. At 



v, 



gw 



< 10 Hz, Vtcwiygw) will continue to fall. Therefore, all constraints fl5.15| ), ( |5.16| ) and 



( p.l?| ) will be satisfied. 



VI. OTHER BACKGROUND NOISES 

The gravitational wave strain amplitude h v due to various galactic binary systems is 
obtained in Ref. KJ (see also |4l|j42]|). Among these binary systems, the background due to 



Close White Dwarf Binaries (CWDBs) is dominant at relevant frequencies, 10 4 Hz < v gw < 



10 *Hz. We obtain the spectral density due to CWDBs, fif D , from Table. 7 in Ref. |40 



where the space density of CWDBs is assumed to take the maximal theoretical value in the 
calculations of Webbink |43[ . 



In addition to this background, there are instrumental noises. The instrumental noises 
can be classified into two kinds, i.e., the optical-path noise and the acceleration noise. The 
optical-path noise, which includes the shot noise, beam pointing instabilities and so on, 
dominates at high frequencies and is estimated as [HI, 



K pt = 4 X x ^ - i2nLUgJc)2 ^ 4 x + L09? (-^-J [Hz-n 

(6.1) 



2L v v gwi ' V VlO" 2 Hz 




where we adopted 2L = 10 10 m. We include the last factor Jl + {2-KLv gw j 'c) 2 considering 
the fact that the sensitivity decreases when the wave period is shorter than the round- 
trip light travel time in one of LISA's arms. The acceleration noise, which includes the 
thermal distortion of spacecraft, the gravity noise due to spacecraft displacement and so on, 



dominates at low frequencies and is estimated as [44 



- 6xlQ" 15 (^/10- 4 Hz)- 1/3 ms-VvTfe ( ^ 1/2 

K ~ 2L(2nu gw y - 1.52 x 10 [Hz j. 

(6.2) 

In Fig. 6, setting M BH = 0.5M© and Vth 2 = 0.1, we summarize the spectral density 
for the cosmological background, h c ° s , in Eq. (|5.9|) and that for the halo background, h^ al °, 
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in Eq. ( |4.15|) together with the spectral density for CWDBs, h]^ D , and that due to the 
two kinds of instrumental noises, /i° pt and h^ cc . When we plot the halo background case, 
we assumed an isothermal halo. As we can see from Fig. 6, the cosmological background 
dominates the halo background for an isothermal halo with this choice of model parameters. 
Since the dependence of the cosmological background on Qh 2 is different from that of the 
halo background, which we can see in Table I and Table II, the halo background dominates 
the cosmological background when the density parameter of BHMACHOs is very small. 
However, since such a situation is less interesting for us, we neglect the halo background 
hereafter. Of course, if the halo is more concentrated, i.e., A is larger and R a is smaller, the 
halo background can dominate the cosmological background, as we can see from Fig. 2. 

At first glance, the spectral density due to CWDBs, h^ D , dominates the cosmologi- 
cal background due to BHMACHOs, h c ° s . However, the cosmological background due to 
BHMACHOs, hl° s , can be more important. This is because the confusion noise level for 
detection of individual sources also depends on the number of noise sources in a frequency 
resolution bin. At lower frequencies many galactic CWDBs exist in a frequency resolution 
bin, while the expected number at higher frequencies becomes less than one. Hence the 
confusion noise level determined by CWDBs decreases suddenly at a transition frequency 
v gw ~ 10~ 2 ' 5 Hz [ |45|j . Therefore, as for the confusion noise level, the cosmological back- 
ground produced by BHMACHO binaries dominates the background due to CWDBs above 
the transition frequency. 



VII. INDIVIDUAL SOURCE 
A. expected number 

Let us consider the expected number of observable individual sources. The frequency 
resolution is determined by the observing time T [f46| , |35| . Here we assume that the frequency 
resolution is given by 

^=f, (7-1) 

with T ~ 1 yr. For the noise amplitude, h u , the threshold amplitude of gravitational waves 
from individual sources is given by, 

h th = 5h u Au 1/2 , (7.2) 

where we set the signal-to-noise ratio (SNR) as 5. [] 

The amplitude of gravitational wave at v gw from a binary with e, p and the distance 
from the earth, d, is given by 



4 Here we do not consider a reduction factor due to the antenna pattern of an interferometer. 
Hence, we overestimate the number of individual sources in the following discussion. A traditional 
way to take this factor into account is to adopt the root-mean-square value of the signal averaged 
over the entire sky, which implies that the reduction factor is \/5 [p5 , 44 1 . 
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c 3 / 2 y^K v gw ' ' 4ird 2 



(7.3) 



where L^(z/ S „,,e) is given by Eq. (O). Then, the requirement that the signal exceeds the 
threshhold, h > h th , determines the maximum distance to individually observable sources, 



D {p) (v e) 



G 



nc 3/2 h th 



(7.4) 



where L$ is given by Eq. ( |4.3| ). Then the expected number of individual sources per unit 
logarithmic frequency in the p-th harmonic is given by 



N (p \i 



gwj 



deJ\f {p) (v gw ,e)v gw -^' 



'gw J P 



and 



n(r) d 3 x. 



(7.5) 



(7.6) 



where j v ^tiy vi e; to) and n(r) are given in Eqs. ( |3.10| ) and ( |4.6| ), respectively. The expected 
number of individual sources per unit logarithmic frequency is given by 



N(u gv) ) = y £N^{v gv> ). 
P =i 

The total expected number in the range of V\ < u gw < u 2 is 



(7.7) 



N 



tot 



lni/2 



In v\ 



N(u gw )d(\nu gw ). 



(7.8) 



When D^ ax (u gw , e) is much smaller than the distance to the center of our galaxy from 
us, D , under the assumption of the uniform density we can approximate Eq. ( |7.6| ) as 



N {p \u, 



47m(A)) 



gu>) 



[D%l(v gw ,e)r = K 



C 2 T LMC 

3GM BH D 2 



(7.9) 



where 



K 



Dl/D 



[i + {di/ddy 



" dy- 



1 + (§- 2 w a vy + y' 



(7.10) 



is a factor obtained from Eqs. ( [4.6|) and (|4.9|) which depends only on the shape of our galaxy. 
In Fig. 7, we plot K as a function of the core radius, R a , in units of kpc for A = 1, 1.5 and 
2. K is about 0.8 when A ~ 1. 
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When such an approximation is not valid, we have to integrate Eq. ([7.61) . When A is 
integer, however, we can integrate Eq. ( |7.6|) analytically. Evaluating the integrals 







drAur n(r) 

/■•DmL(>9it>,e) + A) 

+ / dr2irr n(r) 



1 - 



D2 + r*-(D(±(v gw ,e)y 



2rD 



(7.11) 



for D £L( u 9v» e ) > A), and 



dr2nr 2 n(r) 



Dl + r*-(D(±(v gw ,e)f 



2rDr 



J Do-D\ n ' ax {y gw ,e) 

for D^ x (i/ gw ,e) < D , we find both expressions ([7.1 1|) and ( |7.12|) give the same results, 



(7.12) 



A 



= 1 : A/^W, e) = f[D%Uv gw , e)fn s { [I A (1 + r ) + I A (1 - r )] + \r\ 

+ 3 -^(l-r 2 + r 2 ) In [ ( 1 + r °) 2 + r ' "' 



A = 2 : A^(*V, e) = y [^L(*V, e)] 3 n s { ^r a 4 [J A (1 + r ) + J x (l - r 



(1 -r ) 2 + r 2 



-^ln 

8r 



(l + r ) 2 + r 2 
(1 - r ) 2 + r 2 



>, (7.13) 



(7.14) 



where 



D a 



Dr 



i ,, — — . / n — — , Ia( x ) — — arctan ( — ^ . (7-15) 

Dmax(Vgw,e) Dmax(v gw ,e) r a Vr a / 

When D$ ax {y gw ,e) is much larger than the core radius, D a , and the distance between the 
center of our galaxy and us, D , Eqs. (|7.13|) and ( [7. 14|) can be approximated as 

A = 1 : Ar^(u gwj e) = AnD 2 a D^ ax {u gw , e)n s : D$ ax (u gU]7 e) » D , D a , (7.16) 
A = 2: M {p) {v gw ,e)='K 2 Dln s : D^ ax (u gw , e) > D , D a . (7.17) 



B. the case when the confusion limit is determined by BHMACHOs and 

instrumental noises 

Since neither the abundance of halo CWDBs nor that of cosmological BHMACHOs has 
not been established, we do not know which background determines the confusion limit of 
individual sources. In this subsection, we assume that BHMACHOs dominate the energy 
density of the present universe, and they also determine the amplitude of the gravitational 
wave background. As we have seen in Sec.|VT|, under this assumption it is natural to think 
that the cosmological background dominates the halo background. Then the total noise 
amplitude is determined by 
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K ot = \I{h™y + {hTy + (/ir) 2 , (7. is) 



gwj 



where hl° s , h° pt and h^ cc are given in Eqs. ( |5.9| ), ( |6.1| ) and (|6.2|), respectively. 

In Fig. 8, we plot the expected number of individual sources in each harmonic 
for p = 1,2,3,4 and 5. The solid lines are the expected number for A = 1, D a = 5 kpc, 
Mbh — O.5M and Qh 2 = 0.1. As we can see in Fig. 8, the second harmonic at high 
frequencies is prominent. Therefore, the estimate of the total number does not affected by 
errors in the estimate of the noise level at low frequencies so much, which errors are partly 
caused by truncating the summation at p = 10 in Eq. (|5.7|) . In Fig. 9, we plot the maximum 
distance of observable individual source defined by Dml x (v gw , e vea k) for p = 1,2,3,4 and 5, 
where e vea k is the eccentricity at which the distribution function f Ufi ,t{ v pi e; to) is maximum. 
As we can see in Fig. 9, the maximum distance D^ ax {y gw , e vea k) extends beyond ~ 10 2 
kpc at high frequencies. This means that almost all halo BHMACHO binaries radiating 
gravitational waves at these frequencies can be observed as individual sources. In Fig. 10, 
we plot the total expected number of the individual sources per logarithmic frequency, 



N(v gw ), in Eq. (|7~7|) with the solid line. The dotted line is the expected number using 
Eq. ( |7.9| ) with K = 0.8, which is valid for D^ ax (u gw , e) <C D . The dashed line is the 



expected number using Eq. (|7.16| ), which is valid for D$ ax {v gW) e) ^> D , D a . From this 



figure, we can see that Eqs. ( |7.9| ) and ( |7.16| ) are good approximations for D^ ax (u gw , e) <C D 
and for D$ ax (i> gw ,e) ^> D , respectively. The total expected number, N tot , in Eq. ( |7.8|) in 
the range of 10" 4 Hz < v gw < lO^Hz is about 792 for A = 1, D a = 5 kpc, M BH = O.5M 
and tth 2 = 0.1. 

Now we show N tot is approximately proportional to M £ ^ 2 J 81/ ' 407 for Mbh ~ M Q . Since 
the second harmonic at high frequencies almost completely determines N tot , it is sufficient 
to consider only p = 2 and e€l case. In this case, from Eq. (^^) with Eqs. (|4.3| ), (|5.11|) 



and Q5.12j ), we find the maximum distance Z^^z/^, 0) is proportional to fg^ % M^^ 222 . 
The expected number per unit logarithmic frequency, N(u gw ), has a maximum value when 
^max{ v gu>i 0) ~ Ah as we can see m Fig- 9 and Fig. 10. Therefore, the frequency v gw at 
which N{v gw ) has a maximum value is proportional to M^ 81 ^ 407 . At this frequency, we can 
use the expression for M^(iy gw , 0) given in Eq. (|7.16|) , which turns out to be proportional 
to M^h with the aid of Eq. ( }4.9| ). On the other hand, using Eqs. ( |3.13| ) and ( |3.14j ), we see 
VgwJdef Uje j(v p ,e;to) is approximately proportional to u gw f Vte (u p ;t ) oc ^ /3 Mb]J 0/u1 oc 

M 12 H /m . Thus, we obtain N^(u gw ) oc M" 281 / 407 from Eq. Q. Since the total expected 
number is approximately determined by the maximum value of iV( 2 ) (vgw), the relation N tot oc 
M^^ 81 ' 407 is obtained. Here we note that N tot does not strongly depend on Mbh for Mbh ~ 
M . Q 



5 This dependence cannot hold for extreme values of Mbh- 
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C. the case when the confusion limit is determined by CWDBs and instrumental 

noises 



In this section we assume that gravitational wave background due to CWDBs and the 



instrumental noises determine the confusion limit. As noted at the end of SecjVT], the 
confusion noise level due to CWDBs decreases suddenly at frequency v gw ~ 10~ 2 5 Hz. 
Thus, even if the confusion noise level is determined by CWDBs at low frequencies, the role 
of CWDBs is replaced with BHMACHOs at high frequencies. Therefore, the background 
due to BHMACHOs must be simultaneously taken into account. However, for simplicity, 
we ignore the background due to BHMACHOs in this subsection. Instead, we neglect the 
sudden decrease in the confusion noise level due to CWDBs, and hence the estimate given 
in this subsection will provide a lower bound for the total expected number. 
Now, the total noise amplitude is determined by 

hf = ^{h WD f + (K pt ) 2 + {h^)\ (7.19) 
where h„ D is taken from Table 7 in Ref. [40| and we use the linear interpolation. 



As before, since the expected number in the second harmonic at high frequencies is 
prominent, the total expected number is almost insensitive to the noise uncertainty at low 
frequencies. In Fig. 11, we plot the maximum distance D^> ax (iy gw , e pea k) for p = 1, 2, 3, 4 and 
5. In Fig. 12, we plot the total expected number of the individual sources per unit logarithmic 
frequency, N(u gw ), in Eq. ( [7.7|) . For comparison, with the dotted line, we plot again the 
total expected number per unit logarithmic frequency, N(i/ gw ), when the background due to 
BHMACHOs is dominant. It is the same plot that is presented by the solid line in Fig. 10. 
The total expected number N tot in Eq. ( |7.8|) in the range of 10~ 4 Hz < v gw < 10 _1 Hz is 
about 329. 



VIII. SUMMARY 

In this paper, we have investigated low frequency gravitational waves emitted by BHMA- 
CHO binaries. We have evaluated the gravitational wave background produced by BHMA- 
CHO binaries in the Milky Way halo. We have also evaluated the cosmological gravitational 
wave background produced by the extragalactic BHMACHO binaries, which is larger than 
the halo background when we assume that the density profile of our halo is isothermal. The 
root-mean-square strain amplitude of the gravitational wave background due to CWDBs 
is possibly larger than that due to the extragalactic BHMACHO binaries. However, once 
the number of CWDBs per frequency resolution bin is taken into account, the cosmological 
gravitational wave background due to extragalactic BHMACHO binaries will determine the 
confusion limit for the detection of individual sources at v gw > 10~ 2 ' 5 Hz. We have found 
that the cosmological gravitational wave background produced by BHMACHO binaries is 
well below the existing observational limits obtained from the pulsar timing, the quasar 
proper motions, the big bang nucleosynthesis and so on. We have found that one year ob- 
servation by LISA will be able to reveal at least several hundreds of BHMACHO binaries 
whose gravitational wave amplitudes exceed the confusion limit when we require the SNR 
to be greater than 5 for detection. If the confusion limit is determined by the noises due 
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to extragalactic BHMACHO binaries and instruments, the expected number of observable 
individual BHMACHO binaries will be about eight hundreds. This suggests that LISA will 
be able to clarify whether MACHOs are primordial black holes or not, together with the 
results of LIGO-VIRGO-TAMA-GEO network H^||. 

At frequencies above ~ lO- 2 - 23 (T/lyr)- 6 / 11 (M SH /O.5M )- 5 / 11 , LISA can measure dis- 
tances to BHMACHO binaries since binaries significantly change their orbital frequencies 
through gravitational wave emission within the observation time T [4~I|. Moreover, the 
LISA's angular resolution is roughly 10~ 3 steradians (equivalently, 3 square degrees) at fre- 
quencies v gw ~ 10~ 2 Hz for SNR = 5 |^7J. Therefore, it will be possible to draw a map 
of the distribution of BHMACHO binaries. Recall that the BHMACHO binaries trace the 
mass distribution of the Milky Way halo. Then, if MACHOs are black holes, it may be 
possible to obtain the mass distribution of our galaxy! Note that the expected number of 
BHMACHO binaries is large enough to obtain a precise map, and the maximum distance to 
observable individual sources, ~ 10 2 kpc, is sufficiently large as shown in Fig. 9 and Fig. 11. 
(See Ref. (48| for further discussions.) 

In reality, the mass of MACHOs inferred from observations depends on the halo model 
(e.g. Ref. |13[)- However the total expected number N tot of individually observable BH- 
MACHO binaries does not strongly depend on Mbh fc> r Mbh ~ M Q . In the case of a 
general distribution of MACHOs masses, we should consider binaries made from different 
mass BHMACHOs. This is an interesting future problem. 

There may be other observational signals of BHMACHOs except for gravitational waves 
from BHMACHO binaries which are formed in the early universe. For example, we may 
be able to detect gravitational waves from BHMACHOs orbiting a massive black hole in a 
galactic nuclei [49 -[51], though the event rate is uncertain. Moreover, massive black holes 
in galactic nuclei may be formed from the BHMACHOs through their merging process j52 



Gravitational waves from their mergers may be detected by LISA, though the event rate is 
also uncertain. It is a future problem to obtain a reliable estimate of the event rate for such 
events. 

There are two microlensing events toward the LMC and the Small Magellanic Cloud 
)6l. In both events, the lens objects are not likely in our halo. However, the fraction of 



BHMACHOs that are in binaries with a ~ 2 x 10 14 cm and M BH = 0.5M Q is found to be 
~ 0.26% for Qh 2 = 0.1, taking into account the evolution of the BHMACHO binaries. (In 
Ref. [p7| , ^| ) we did not take the evolution into account). Therefore, no detection of halo 
binary lens events does not conflict with the BHMACHO scenario. 
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APPENDIX A: THE CONTRIBUTION FROM HIGH P-TH HARMONICS 



In this section, we investigate how large the contribution from high p-th harmonics to 
the gravitational wave flux in Eq. (f4.14|) . As we can see in Eq. ( |4.14| ), the gravitational wave 
flux from BHMACHO binaries in the Milky Way halo per unit frequency at v gw is given by 

oo 

FuM = il^Y.Hp), (ai) 
p =i 

where 

F(p) = P~ 13/3 I de f VA t{u Pi e; t )g(p, e). (A2) 



For a fixed large p > 100, g(p,e) has a peak at e ~ 1, and the maximum of g(p,e) is 
about several times of g(p, 1). f\ Therefore, for a fixed large p ^ 100, g(p,e) has an upper 
bound as 



p 2 p 2 

gip, e) ^ g(p, l) = -^Jpip) ~ g- 



r(i/3) - 



22/331/6^1/3 



(A3) 

27/334/3^2^ ' ^ 



where we use Eq. (Al) in Ref. |32j at the second equality and we use the asymptotic formula 
for Bessel function at the third equality. Using this relation, Eq. ( |A2| ) becomes 



T{p)=p 13/3 / def u , e)t (v p ,e;t )g( y p,e)<p 13/3 g(p,l) f de f v , e ,t(v P , e; t ). (A4) 
J Jo 

Next, we integrate 

/ def u e j(v p ,e;t ) = — de f(a p ,e;t ), (A5) 

JO 6U p JO 

in Eq. (|A4j) where the orbital frequency v v = v gw jp is given in Eq. and the semimajor 
axis a p is given in Eq. (|3.9|). [] The distribution function at the present epoch, f(a, e; t ), is 
given by Eq. (|3.S| ) and the initial distribution function, /j(aj,ej), is given by Eq. (|2if ). The 
ranges of ai and are < < 00 and -^1 — P 2 < < 0, respectively. 

Since e) has a peak at e ~ 1 for a large p, it is sufficient to consider only 1 - e 2 < 1 
case. When eccentricity is near 1, we can express and by a p and e as 

Oj = a p w~ 2 [l + u} 2 , (A6) 
l-e 2 = (l-e 2 ) M 2 [l + u]- 2 , (A7) 

where 



6 We have not proved this fact. However we confirmed this numerically. 

7 We add the subscript p to a in order to emphasis that the semimajor axis a depends on p for a 
fixed gravitational wave frequency v gw . 
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u 



a 



) 4 (l-e 2 ) 7 / 2 . 



(A8) 



!9|) and the fact that the difference 



The above equations can be derived from Eqs. 
between the coalescence time at t = tf and that at t = to is to. For given a p , there is e p that 
satisfies u = 1 as 



-8/7 



(A9) 



For e p < e < 1, we can approximate Eqs. ( |A6[ ) and ( |A7| ) as 



-2 



(A10) 



since u < 1. For < e < e„, we can approximate Eqs. 



and ( |A7| ) as 



e 2 ~ 



(All) 



since u > 1. Note that Eq. ( All ) is a good approximation for w 3> 1 even when 1 — e 2 ~ 1. 

We separate the integration range of Eq. ( |A5|) into two parts, < e < e p and e p < e < 1. 
In each range, we use the approximate relations of Eqs. ( |A10|) and (|A11|) . Firstly, we 



integrate Eq. 



in the range of < e < e p . Since e ~ ej for < e < e p and VI — /3 2 < 



ej < 1 have to be satisfied, the integration range is y/1 — (3 2 < e < e p . The present 
distribution function f(a p , e;t ) is the same as the initial distribution function /^(a^e^), 
since the relations in Eq. ( |A11| ) are satisfied for < e < e p . We perform the integration as 



/ P de fv,e,t{v P , e; t Q ) = f " de f{a p , e; t 
Jo iv v Jo 



2 a p 
3z/„ 



2v p 
1 

2v n 



V 



ax J 



def(a p , e;t ) ~ ^ 



defi(a. 







'_ 






n 


exp 




— exp 






"_ f^E. " 


f 






exp 




— exp 



/->■ 







a p \^ 



fOp\ 

\a J 



4/7 



(1 -e2)V2 

a 



ax 



3/4 




J_ 

2;/, 
2 





-4/7, 



\ax ) \aoJ 
u p \axj 



a p < (3 ' ao 
/3- 4/7 ao < a P £ /3- 28 / 37 aJ 6/37 (ax) 21 /37 

/5 _ 28/3 7 a 16/37 (a - )21/37 < ^ < a - 

ax < a„ 



2^te) 3/2 te) 4/ v 

Ivgw \ ax J r 



4/7 



50/21 



3/2 



P < Pi 

Pi < P ^ P2 

P2^P<P3 
P3 < P, 



(A12) 



where a gw is defined as 
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"gw 



(27TI/j 



gw) 



-2/3 



[G(M BH i + M BH2 )} 



1/3 



P 



-2/3 



"i 



I" 



(A13) 



and 



Pi 

P2 
P3 



p-6/7 I a gw 

^-42/37 f a gw \ 24//37 / O'gw 



-3/2 



Qgw 

ax 



a J 

-3/2 



-63/74 



(A14) 
(A15) 
(A16) 



In Eq. ( |A12[ ), the first range of p is determined by the condition, e p < a/1 — /3 2 , where 
the integration range becomes zero. The second range of p is determined by the condition, 
P(a p / ao) 4//7 {dpi ax) 3//4 < 1, where we can expand out the second exponential in Eq. QA12 ) 
with good accuracy. The third range of p is determined by the condition, (a p /ax) 3 ^ 4 = 1, 
where the integration begins to damp exponentially. 

Secondly, we integrate Eq. ( |A5| ) in the range of e p < e < 1. Using the relations in 
Eq. ( |A10| ), we perform the integration as 



2a r , 



3u„ 



2a r 



def v ^t(v p ,e]t Q ) = — / def(a p ,e;t ) ~ — / de — 



3z/„ 



1 



,2\ 3/2 



1 , 1 



1/2 . /i _ 



l-e 2 ) 

2\ 3/2 



x exp 
"i 



2^ p (l-e 2 ) 3 / 2 [axfl 2 a p \l - e 2 
(3 fa,;\m 



ede 



{l-e 2 ) 1 / 2 \ax 
1 (3 



2v p {\ 



ede ■ 



3/2 



3/2 



ax J 



u exp 







2U/2 



e 2 ) 



ax 



3/4 
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-5/2 
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-12 



x exp 
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37za, 
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-4/37 
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37z/ 
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45/74 
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-4/37 / 44x 
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4/3 



Using Eqs. ( [All ) and flATTD , .F(p) in Eq. 



\37) P 
can be approximated as 



(A17) 



Hp) £ P~ 13/3 g(p, i) /&/, 



F(l/3)] z 3 

2 7/3 3 4/ 37r 2^ 



de f u>e ,t(vp, e; t ) + / de L e Av p , e; £ ) 
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(A18) 



The summation of ^(p) with respect to p in Eq. (|A1|) can be approximated as the integration 
with respect to p. The integration of J^{p) can be performed as 
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+ fW (A) 8 "" ( ± _ (A) 3 ' 2 _ ^ (4) 3/2 (^fpH, (A19) 
lb Vao;/ Vao/ Vao;/ lb Vax/ Vao/ I 



where we consider the case of pi < Po < P2- All the contribution from high p-th harmonics 
with po < p to the gravitational wave flux in Eq. (|AT|) is given by 



,,„., < iLy g °j 3 s( Po ). 



In Fig. 13, we plot F^\u gw ) / F^\u gw ) for M BH = O.5M , Mi 2 
where F^ po \v gw ) is the gravitational wave flux summed up till p 
p(p<po) (y gw } is less than a few percent at v gw < 10~ 25 Hz. At v gw > 10~ 2 ' 5 Hz, the error 



(A20) 

0.1 and po = 1000, 
= po- The error of 



of Fjv^\p gw ) seems to be large, however this is merely overestimation of the error. As we 
can see in Fig. 3, the results summed up till p = 10 coincide with the results summed up till 
p = 1000 very well at v gw > 10~ 2,5 Hz, and this indicates that the error due to the cutoff 
of high harmonics is small at v gw > 10~ 2,5 Hz. To conclude, it is sufficient to sum up till 
p ~ 1000 in Eq. ( |Al| ) for a few percent precision. 
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FIGURE CAPTION 



Fig.l. The probability distribution function f Vfi A v vi e > *0 as a function of eccentricity e for 
(a) v v = 10" 1 Hz, (b) u p = KT 2 Hz, (c) v v = 10~ 3 Hz and (d) v v = 10" 4 Hz, at present 
epoch, i.e. we set t 1 = t Q = 10 10 yr. Here we adopt M BH = O.5M and flh 2 = 0.1. 



Fig. 2. The shape factor I defined in Eq. ( |4.13| ) as a function of the core radius D a for three 
values of A. 

Fig. 3. The solid line is the spectral density of the gravitational wave strain amplitude 
fohaioj^fj f Qr Yiq^q background gravitational waves for M BB = O.5M and VLh 2 = 
0.1. The long dashed line is the spectral density obtained by the summation up to 
p = 10 in Eq. ( 4.14|) . The dotted line is the approximate spectral density, replacing the 



flux F v in Eq. ( |4.15| ) with the approximate flux F v in Eq. ( |4.1(j| ). The approximation 
is valid for e < 1. Note that h^ alo /\I does not depend on the halo shape factor, /. 



Fig. 4. The solid line is dVL GW / 'd(ln R) in Eq. (|5.8|) as a function of R/Rq for M BH = O.5M 
and VL BH h 2 = 0.1 at four frequencies. The dotted line is the approximation of 
dft GW I d(\n R) , replacing the luminosity C G w i n Eq. ( |5.7|) with the approximate lumi- 
nosity Cqw m Eq. ( |5.10| ). The approximation is valid for e 1. Solid and doted lines 
are almost same for high frequencies. 

Fig. 5. The energy density of the cosmological gravitational wave background per logarithmic 
frequency £lGw{v gw ) as a function of the gravitational wave frequency v gw , where we 
summed up to p = 10 in Eq. (|577| ) . We adopt M B h = O.5M and VLh 2 = 0.1. 

Fig. 6. The spectral density h v for the cosmological gravitational wave background in 
Eq. ( |5.9|) by the solid line, for the halo background in Eq. ( |4.15| ) by the dashed line, 
the CWDBs background by the cross, and the instrumental noise in Eqs. ( |6.1| ) and 
( |OD by the dotted line. We adopt M BH = O.5M and Qh 2 = 0.1. The shape factor / 
in Eq. ( [4.13[ ) for the halo background is taken to be 7. 

Fig. 7. The shape factor K defined in Eq. ( |7.10|) as a function of the core radius D a for three 
values of A. 

Fig. 8. The expected number of the individual sources per logarithmic frequency in each 
harmonic N^ p \u gw ) in Eq. (|7.5| ) for p = 1,2,3,4,5, We adopt A = 1, D a — 5 kpc, 
M BH = O.5M and VLh 2 = 0.1. 

Fig. 9. The maximum distance D^l x (iy gw , e pea k) to a source whose gravitational wave am- 
plitude exceeds the threshold value for p = 1,2,3,4,5, where e pea k is the eccentricity 
at which the distribution function f v ^A v vi e ' *o) is maximum. The background is due 
to cosmological BHMACHO binaries and instruments. We adopt M BB = O.5M , 
VLh 2 = 0.1, A = 1 and D a = 5 kpc. 

Fig. 10. The total expected number of the individual sources per logarithmic frequency 
N^gw) in Eq. ( |7.7|) . The background is due to cosmological BHMACHO binaries 
and instruments. The dotted line is the expected number using Eq. (|7.9| ) with K=0.8, 
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which is valid for D$ ax (i> gw ,e) <C D . The dashed line is the expected number using 
Eq. ( fnq ), which is valid for D$ ax (v gw , e) > D . We adopt M BH = 0.5M Q , Vth 2 = 0.1, 
A = 1 and D a = 5 kpc. 

Fig. 11. The maximum distance D^ ax {y gw ^e vea k) to a source whose gravitational wave am- 
plitude exceeds the threshold value for p = 1,2,3,4,5, where e pea k is the eccentricity 
at which the distribution function fu,e,t( u pj e '^o) is maximum. The background is due 
to CWDBs and instruments. We adopt M BH = O.5M , Qh 2 = 0.1, A = 1 and D a = 5 
kpc. 

Fig. 12. The total expected number of the individual sources per logarithmic frequency 
N{v gw ) in Eq. (fT7| ). Solid line is the case that the background is due to CWDBs 
and instruments, and dotted line is the case that the background is due to cosmologi- 
cal BHMACHO binaries and instruments. We adopt M BH = O.5M , Vth 2 = 0.1, A = 1 
and D a = 5 kpc. 

Fig.13. F^\u gw )/F^\u gw ) for M BH = O.5M , VLh 2 = 0.1 and p = 1000, 
where F^°\v gw ) is the gravitational wave flux summed up till p = po in 
Eq. ( [4.14j ) or Eq. (|A1|) . The dotted line is the three percent error line, i.e., 



F^\v gw )/F^\v gw ) = 0.03. 
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TABLES 



TABLE I. The spectral density of the gravitational wave strain amplitude h^ al ° /I 1 / 2 [Hz l / 2 } 



for the halo gravitational wave background in Eq. (4.15). The shape factor I is shown in Fig. 2 



i „ / n i l \ 
log(z/ slu [HzJ) 






i ™ ( l halo 1 t\ /2 

log(n,J; /i ' 


[Hz / j) 








M B h = 0.5M Q 




7\ /f a at j f 

Mm = 0.05M t 


7) 


M B H = 5M 






\ltl = 0.1 


\ln = 0.5 


n[,2 n 1 
SZ/l = 0.1 


O i 2 A r 


\lh = 0.1 


iZ/l = 0.5 


-1.00 


-21.18 


-21.03 


-21.58 


-21.43 


-20.78 


-20.63 


1 OA 

-1.20 


-20.95 


OA OA 

-20.80 


o i or 

-21.35 


1 OA 

-21.20 


oa r r 

-20.55 


OA 1 A 

-20.40 


-1.40 


OA r?ri 

-20.72 


oa r r* 

-20.56 


-21.11 


oa a/ 1 

-20.96 


OA O 1 

-20.31 


oa 1 r* 

-20.16 


-1.60 


-20.48 


-20.33 


-20.86 


-20.72 


-20.08 


-19.93 


1 on 

-1.80 


oa o r 

-20.25 


OA 1 A 

-20.10 


OA C O 

-20.58 


OA A f* 

-20.46 


1 a or 

-19.85 


1 A 

-19.70 


-2.00 


-20.01 


-19.86 


-20.32 


-20.18 


-19.61 


-19.46 


-2.20 


1 A 

-19.77 


1 A /' ("1 

-19.62 


OA 1 ^ 

-20.14 


1 A A A 

-19.94 


1 A O O 

-19.38 


1 A OO 

-19.23 


-2.40 


1 a r a 

-19.50 


1 A OT 

-19.37 


OA AO 

-20.08 


1 A OO 

-19.82 


i a i r 

-19.15 


1 A AA 

-19.00 


-2.60 


1 A OO 

-19.22 


-19.09 


-20.10 


-19.79 


-18.91 


-18.76 


-2.80 


-19.01 


-18.83 


-20.17 


-19.83 


-18.67 


-18.52 


-3.00 


-18.90 


-18.66 


-20.23 


-19.90 


-18.42 


-18.28 


-3.20 


-18.89 


-18.60 


-20.25 


-19.95 


-18.14 


-18.01 


-3.40 


-18.95 


-18.62 


-20.23 


-19.95 


-17.89 


-17.73 


-3.60 


-19.02 


-18.68 


-20.20 


-19.93 


-17.74 


-17.52 


-3.80 


-19.06 


-18.74 


-20.17 


-19.90 


-17.70 


-17.42 


-4.00 


-19.05 


-18.76 


-20.15 


-19.87 


-17.73 


-17.41 


-4.20 


-19.02 


-18.74 


-20.14 


-19.85 


-17.80 


-17.46 


-4.40 


-18.99 


-18.71 


-20.13 


-19.84 


-17.85 


-17.53 



32 



TABLE II. The spectral density of the gravitational wave strain amplitude h c ° s [Hz 1 I 2 ] for 
the cosmological gravitational wave background in Eq. (5.9). 

log(^[Hz]) log(7,;."*[Hz '/-]) 
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